home *** CD-ROM | disk | FTP | other *** search
- -- -----------------------------------------------------------------------
- -- Title: primitive_functions
- -- Last Mod: Thu Apr 19 11:25:43 1990
- -- Author: Vincent Broman <broman@nosc.mil>
- -- Copyright 1990 Vincent Broman
- -- Permission granted to copy, modify, or compile this software for
- -- one's own use, provided that this copyright notice is preserved intact.
- -- Permission granted to distribute compiled binary copies of this
- -- software which are linked in with some other application.
- -- Permission granted to distribute other copies of this software,
- -- provided that (1) any copy which is not source code, i.e. not in the
- -- form in which the software is usually maintained, must be accompanied
- -- by a copy of the source code from which it was compiled, and (2) the
- -- one distributing it must refrain from imposing on the recipient
- -- further restrictions on the distribution of this software.
- --
- -- Visibility:
- -- Description:
- -- functions on floating point types whose implementation
- -- depends on access to the mantissa/exponent representation
- -- of the floating point number. this includes
- -- integer/fraction operations.
- --
- -- This is a portable, slow implementation.
- --
- -- Exceptions: numeric_error upon overflow of function scale,
- -- constraint_error only if type Float is constrained.
- -- -----------------------------------------------------------------------
- package body primitive_functions is
- --
-
-
- subtype expbits is integer range 0 .. 6; -- 2**(2**7) might overflow
-
- Float_exp: array( expbits) of integer :=
- (1, 2, 4, 8, 16, 32, 64);
- Float_power: array( expbits) of Float :=
- (2.0 ** 1, 2.0 ** 2, 2.0 ** 4, 2.0 ** 8,
- 2.0 ** 16, 2.0 ** 32, 2.0 ** 64);
- Float_neg_power: array( expbits) of Float :=
- (0.5 ** 1, 0.5 ** 2, 0.5 ** 4, 0.5 ** 8,
- 0.5 ** 16, 0.5 ** 32, 0.5 ** 64);
-
- -- double stuff deleted here
-
- function exponent( x: Float) return integer is
- --
- -- return the exponent k such that 1/2 <= x/(2**k) < 1,
- -- or zero for x = 0.0 .
- --
- begin
- if x = 0.0 then
- return 0;
- end if;
- -- nonzero x
- -- essentially, just multiply by powers of two to get in range
-
- declare
- ax: Float := abs( x);
- exp: integer := 0;
- -- ax*2.0**exp is invariant
- begin
- while ax >= Float_power( expbits'last) loop
- ax := ax * Float_neg_power( expbits'last);
- exp := exp + Float_exp( expbits'last);
- end loop;
- -- ax < 2**64
- for n in reverse expbits'first .. expbits'last - 1 loop
- if ax >= Float_power( n) then
- ax := ax * Float_neg_power( n);
- exp := exp + Float_exp( n);
- end if;
- -- ax < Float_power( n)
- end loop;
- -- ax < 2
- while ax < Float_neg_power( expbits'last) loop
- ax := ax * Float_power( expbits'last);
- exp := exp - Float_exp( expbits'last);
- end loop;
- -- 2**-64 <= ax < 2
- for n in reverse expbits'first .. expbits'last - 1 loop
- if ax < Float_neg_power( n) then
- ax := ax * Float_power( n);
- exp := exp - Float_exp( n);
- end if;
- -- Float_neg_power( n) <= ax < 2
- end loop;
- -- 1/2 <= ax < 2
- if ax < 1.0 then
- return exp;
- else
- return exp + 1;
- end if;
- end;
-
- end exponent;
-
-
- function mantissa (x: Float) return Float is
- --
- -- return scale( x, - exponent( x)) if x is nonzero, 0.0 otherwise.
- --
- begin
- if x = 0.0 then
- return 0.0;
- end if;
- -- nonzero x
- -- essentially, just multiply by powers of two to get in range
-
- declare
- ax: Float := abs( x);
- begin
- while ax >= Float_power( expbits'last) loop
- ax := ax * Float_neg_power( expbits'last);
- end loop;
- for n in reverse expbits'first .. expbits'last - 1 loop
- if ax >= Float_power( n) then
- ax := ax * Float_neg_power( n);
- end if;
- end loop;
- while ax < Float_neg_power( expbits'last) loop
- ax := ax * Float_power( expbits'last);
- end loop;
- for n in reverse expbits'first .. expbits'last - 1 loop
- if ax < Float_neg_power( n) then
- ax := ax * Float_power( n);
- end if;
- end loop;
- if ax >= 1.0 then
- ax := ax * 0.5; -- 1 <= ax < 2
- end if;
- if x > 0.0 then
- return Float( ax);
- else
- return Float( - ax);
- end if;
- end;
-
- end mantissa;
-
-
- function scale (x: Float;
- k: integer) return Float is
- --
- -- return x * 2**k quickly, or quietly underflow to zero,
- -- or raise an exception on overflow.
- --
- begin
- if x = 0.0 or k = 0 then
- return x;
- end if;
- -- nonzero x and k.
- -- essentially, just multiply repeatedly by 2**(+-2**n).
- -- if is_
- declare
- y: Float := x;
- exp: integer := k;
- -- y*2.0**exp is invariant
- begin
- while exp <= - Float_exp( expbits'last) loop
- y := y * Float_neg_power( expbits'last);
- exp := exp + Float_exp( expbits'last);
- end loop;
- for n in reverse expbits'first .. expbits'last - 1 loop
- if exp <= - Float_exp( n) then
- y := y * Float_neg_power( n);
- exp := exp + Float_exp( n);
- end if;
- end loop;
- -- exp >= 0
- while exp >= Float_exp( expbits'last) loop
- y := y * Float_power( expbits'last);
- exp := exp - Float_exp( expbits'last);
- end loop;
- for n in reverse expbits'first .. expbits'last - 1 loop
- if exp >= Float_exp( n) then
- y := y * Float_power( n);
- exp := exp - Float_exp( n);
- end if;
- end loop;
- -- exp = 0
- return Float( y);
- end;
-
- end scale;
-
- --
- -- truncation requires making all fractional bits zero in a signed
- -- magnitude representation of floating point numbers.
- -- for large numbers, this has no effect, as all large floats are integers.
- -- for abs( x) < 1.0 this produces zero.
- --
-
- function truncate (x: Float) return Float is
- --
- -- truncate x to the nearest integer value with absolute value
- -- not exceeding abs( x). No conversion to an integer type
- -- is expected, so truncate cannot overflow for large arguments.
- --
- --
- large: Float := 1073741824.0;
- type long is range - 1073741824 .. 1073741824;
- -- 2**30 is longer than any single-precision mantissa
- rd: Float;
-
- begin
- if abs( x) >= large then
- return x;
- else
- rd := Float ( long( x));
- if x >= 0.0 then
- if rd <= x then
- return rd;
- else
- return rd - 1.0;
- end if;
- else
- if rd >= x then
- return rd;
- else
- return rd + 1.0;
- end if;
- end if;
- end if;
- end truncate;
-
-
- function shorten (x: Float;
- k: integer) return Float is
- --
- -- set all but the k most significant bits in the mantissa of x to zero,
- -- i.e. reduce the precision to k bits, truncating, not rounding.
- -- shorten( x, k) = 0.0 if k < 1 and
- -- shorten( x, k) = x if k >= Float'machine_mantissa.
- --
- kex: constant integer := k - exponent( x);
-
- begin -- shorten
- -- the tests on k are needed only to avoid under/overflow
- if x = 0.0 or k < 1 then
- return 0.0;
- elsif k >= Float'machine_mantissa then
- return x;
- else
- return scale( truncate( scale( x, kex)), - kex);
- end if;
- end shorten;
-
-
- function odd (x: Float) return boolean is
- --
- -- predicate indicates whether or not truncate( x) is an odd integer.
- --
- x2: constant Float := truncate(x) / 2.0;
-
- begin
- return x2 /= truncate( x2);
- end odd;
-
- --
- -- all the other integer/fraction functions are based on truncate.
- --
-
- function floor (x: Float) return Float is
- --
- -- return as a Float the greatest integer value <= x.
- --
- sx: constant Float := x;
- tx: constant Float := truncate( sx); -- no constraint_error
- begin
- if sx >= 0.0 or sx = tx then
- return Float( tx);
- else
- return Float( tx - 1.0); -- exactly
- end if;
-
-
- end floor;
-
-
- function ceiling (x: Float) return Float is
- --
- -- return as a Float the least integer value >= x.
- --
- sx: constant Float := x;
- tx: constant Float := truncate( sx); -- no constraint_error
- begin
- if sx <= 0.0 or sx = tx then
- return Float( tx);
- else
- return Float( tx + 1.0); -- exactly
- end if;
-
-
- end ceiling;
-
-
- function round (x: Float) return Float is
- --
- -- return as a Float the integer value nearest x.
- -- in case of a tie, prefer the even value.
- --
-
-
- sx: Float := x;
- tx: Float := truncate( sx);
- diff: Float := sx - tx; -- computed exactly
- begin
- if abs( diff) < 0.5 then
- return Float( tx);
- elsif diff > 0.5 then
- return Float( tx + 1.0);
- elsif diff < -0.5 then
- return Float( tx - 1.0);
- elsif diff = 0.5 then
- if odd( tx) then
- return Float( tx + 1.0);
- else
- return Float( tx);
- end if;
- else -- diff = -0.5
- if odd( tx) then
- return Float( tx - 1.0);
- else
- return Float( tx);
- end if;
- end if;
-
-
- end round;
-
- end primitive_functions;
-
- -- $Header: p_primitive_functions_b.ada,v 3.13 90/04/19 12:34:22 broman Rel $
-